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Abstract Dual Scattering Channel schemes generalise Johns' TLM algorithm 
and replace the latter in situations where the transmission line picture of wave 
propagation fails. This is notoriously the case in applications to fluid dynamics, for 
instance. In this paper, a DSC numerical solution of the Oberbeck-Boussinesq 
equations is presented, which approximate the Navier-Stokes equations for vis- 
cous quasi incompressible flow with moderate variation in temperature. 
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1. Introduction 

Crushing the continuum down into mesh cells is a queer, artificial exercise. 
There are yet natural ways of computing the fields in a cellular mesh, and so 
to mitigate the desastcr. Imposing a cellular mesh is tantamout to locally 
enforcing cell-boundary duality upon space - and the DSC setup offers a 
natural way to handle such situations. 

Dual Scattering Channel (DSC) schemes are characterized by a two- 
step cycle of iteration which alternately updates the computed fields within 
cells and on their interfaces. If the updating instructions arc explicit, then 
a near-field interaction principle leads to the typical structure of the DSC 
algorithm, viz. to its scattering process interpretation. A pair of vectors 
which represent the same field within a cell and on its surface essentially 
constitutes a scattering channel. An equivalent definition can be given in 
terms of a pair of distributions that 'measure' the field within the cell and 
on one of its faces. In the primal DSC scheme which is the Transmission 
Line Matrix (TLM) method along Johns' line |JoB| ports of transmission 
line links visualise these distributions (finite integrals, in this case). 

DSC schemes and their relations to the TLM method [Ch, Tlml-3, Re] 
have been conceptually and technically scrutinised in-depth |Hel) . They are 
unconditionally stable under quite general circumstances - made tangible 
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with the notion of a-passivity |He2) . and they are especially suitable for 
handling boundary conditions, non-orthogonal mesh, or also for replacing a 
staggered grid where otherwise need of such is. 

In section we first recapitulate some characteristic features of DSC 
schemes, which in extenso are treated in |Hel| . before we deal in sectional 
with the Oberbeck - BOUSSINESQ approximation to the Navier-Stokes 
equations. The DSC model outlined here should be considered a prototype, 
first of all. In fact, the BOUSSINESQ equations for viscous fluid flow, in- 
spitc of retaining the usually predominant non-linear advective part of the 
Navier-Stokes momentum equations, can only claim limited range of va- 
lidity, due to their well known simplifications. The Oberbeck-Boussinesq 
approach is, however, prototypical also in providing a basis for many turbu- 
lence models |ATP| that can be implemented following essentially the lines 
of this paper. 

The treated implementation with unstructured hexahcdral mesh, out- 
lined in section^ has recently been implemented and coupled to Spinner's 
Maxwell field solver, thus allowing for computing conductive and convective 
heat transfer simultaneously with the electric and magnetic heat sources in 
a lossy Maxwell field. Besides the underlying ideas and simplifications that 
enter the Oberbeck-Boussinesq approximation and its DSC formulation, 
some numerical results are displayed to illustrate and validate the approach. 



2. Elements of DSC schemes 

In this section, we resume some typical features of DSC schemes, referring 
the technically interested reader always to the systematic exposition |Hel| . 

Given a mesh cell, we think of a port as a vector valued distribution, 
associated to a cell face (with support, however, not necessarily confined 
to that face, cf. scct0J, which assigns a state vector z p = (p, Z) to a 
physical field Z , the latter represented by a smooth vector valued function 
in space-time. We also require that in the given cell a nodal image p~ of p 
exists, such that 

(1) (p~, Z) = (po a, Z) = (p, Z o a- 1 ) , 

for every Z (of class C°° , e.g.) , where a denotes the spatial transla- 
tion a : M 3 — > M 3 that shifts the geometrical node (centre of cell) onto the 
(centre of) the respective face, cf. Figl. 

DSC fields split thus into port and node components, z p and z n , which 
represent the field at the cell surfaces and within the cells. The two compo- 
nents are updated at even and odd integer multiples, respectively, of half a 
timcstep r and are usually as step functions constantly continued over the 
subsequent time intervals of length r . Moreover, we assume that the updat- 
ing instructions are explicit, i.e., with possibly time dependent functions F 
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node V 




and G , for t = mr ; m 6 N 




and we agree upon fixing z p - n (t) := for t < 0. (The 'back in time 
running' form of the sequence has certain technical advantages, cf. [Hell . 1 ) 

A fundamental DSC principle is near-field interaction, which spells that 
every updated state depends only on states (up to present time t ) in the 
immediate neighbourhood. More precisely: The next nodal state depends 
only on states (with their history) in the same cell and on its boundary, and 
a subsequent port state depends only on states (again with history) on the 
same face and in the adjacent nodes. 

As a consequence of near-field interaction, every DSC process allows for 
an interpretation as a multiple scattering process in the following sense. 

If M is a mesh cell system and <9£ denotes the boundary of cell ( £ M, 
then every DSC state obviously permits a unique scattering channel repre- 
sentation in the space 



with canonical projections n^ n : P — > p pn into port and node compo- 
nents of cell £ ( the cell index is omitted, in general, if there is no danger of 
confusion). Also, there is a natural involutary isomorphism nb : P — > P 



which is named the node-boundary map and obviously maps P p onto P n 
and vice versa. For every DSC process z = (z p , z n )(t) , the following 
incident and outgoing fields zf n and z™ ut are then recursively well (viz. 
uniquely) defined, and arc processes in P p and P n , respectively: 



For t < , zf n (t ) : = z% ut (t - § ) : = , and for every < t = mr ; 



P 



nb : (z p , z p ) 




)• 



m e N 



(3) 



z? n (t) := z p {t)-nboz: ut (t- T -), 
*Zt(t+l) ■= z n (t+ \) -nbozl(t). 
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Hence, at every instant holds z p ( t ) = nb o z™ ut (t — J ) + zf n (t) and 
z n (t + ^) = nbo zf n (t) + Zg Ut {t + 5) . Then near- field interaction 
implies that every state is only a function of states incident (up to present 
time t ) on scattering channels connected to the respective node or face. 
Precisely, it is shown that 

Theorem 1 . A pair of functions 31 and C exists, such that for every cell 
(eM the process z" = ir ™ o z complies with 

(4) z?(t + 1) = - /"")) peahen) 
and i/ie port process z £ = 7r ^ o 2 satisfies 

4(t + r) = e((z^(t + I- /iT )) n|af;/ieN ) . 

(5) 2 

( ' I ' short-hand for 'adjacent to ' ) 

(i) The statements imply, of course, that z" out and Z£ in are themselves 
functions of states incident on connected scattering channels, since 

7- 

z c"°ut( <+ o) = ^((^LC* -MT))pe9C;MeN) - z? in {t) and 

(6) ~ 

= e((s£. t (f - - -Mr)) n | p;/ieN ) - z c P out (t - -) 

(ii) 3? and C are named the reflection and connection maps, respectively, 
of the DSC algorithm. 

(iii) Near field interaction implies computational stability, if the reflection 
and connection maps are contractive or a-passive |He2j . in addition. 

Proof. Merc translation of the near-field interaction principle and induction: 
The statement is trivial for t < . If, for instance, (0J holds until t — Z , 
then in virtue of J5J and remark (i) also 

z n C (t + 5) = F ([ z\ } t , [ z? ] t _ f ) 

*f lin (*) + n6« c % ut (t-5 ), K[*c, in ]t-T 

is a function of states incident from p € 9( until f . 



3. The dynamic equations 

DSC algorithms arc thus simply characterized as two-step explicit schemes 
that alternately update states in ports and nodes of a cellular mesh and 
which, in virtue of a near-field interaction principle, allow for a canonical 
interpretation as scattering processes. The latter exchange incident and re- 
flected quantities between cells and their interfaces. 

The ports and nodes are related to physical fields by vector valued dis- 
tributions that evaluate the fields at cell faces and within the cells of a 



DSC NUMERICAL SOLUTION OF OBERBECK-BOUSSINESQ EQUATIONS 5 



cellular mesh. Such a distribution may be a finite integral, as in the case 
of the TLM method, where finite path integrals over electric and magnetic 
fields are evaluated in a discrete approximation to Maxwell's integral equa- 
tions |He3| . In the simplest case, it is only a Dirac measure that pointwise 
evaluates a field ( or a field component ) within a cell and on its surface. The 
distribution can also be a composite of Dirac measures that evaluate a field 
at different points in the cell - which applies, for instance, to the gradient 
functional treated in sect.0J 

Classical thermodynamics with, in particular, energy conservation entail 
the convection- diffusion equation for the temperature T in a fluid of veloc- 
ity u with constant thermal diffusivity a , heat sourcc(s) q , and negligible 
viscous heat dissipation, viz. 

dT 

(7) — h u grad T = a AT + q . 

at 

This is the energy equation for Boussinesq-incompressible fluids, e.g. |GDN| . 
The Navicr-Stokes momentum equations for a fluid of dynamic viscosity fi , 
under pressure p , and in a gravitational field of acceleration g require 

d 

(8) — ( £>u) + ( u • grad )(gu) + gradp = fj, A u + q g . 

The Oberbeck-Boussinesq approximation |Obb| . [Ess] starts from the as- 
sumption that the fluid properties are constant, except fluid density, which 
only in the gravitational term varies linearly with temperature; and that 
viscous dissipation can be neglected. Equations (JSJ) become so with goo = 
const and q(T) =g 00 /3(T(t,x) - T x ) ; /3 := g- 1 d g / dT 

(9) |H + (u.H)u+^ = ^Au^(T(i,x)-r m )g. 

at goo goo 

With the Gauss-Ostrogradski theorem applied to the integrals over A = 
div grad on cell C with boundary <9£, equations |Q 0) yield, with a time 
increment r , the following updating instructions for T and u averaged over 
the cell volume Vq 

T{t+\) ■ = 

( 10 ) a f If 
T + t { u ■ gradT H / gradT ■ dF H qdV } 



and 



(11) 



, t . . . ,, qradp , 

u (t+7;) : = u - t { ( u • grad ) u + } + 

2 g a 



T ^TF L - i gradu-dF + 13 (T - Too) g} 

V^Qoo Joe 



At the right-hand sides enter, of course, the last former updates (at time 
t — r/2 and t, respectively) of the nodal and cell face quantities. 
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T and u arc so updated at the reflection step of the DSC algorithm. 
In contrast, the cell surface integrals at the right-hand sides, in particular 
the gradients that enter these, are updated on the connection step. 
The next section proceeds with that in unstructured hexahedral mesh. 



4. The non-orthogonal hexahedral cell 



The physical interpretation of a DSC algorithm associates a smoothly vary- 
ing, i.e. in time and space sufficiently often continuously differentiable ( for 
instance, C°° ) scalar or vector field Z to port and node states z p and 
z n of a mesh cell system. 

For notational economy ( so avoiding many ' s ) in the following we 
adopt Einstein's convention to sum up over identical right-hand sub and 
superscripts within terms where such are present, while summation is not 
carried out whenever a sub or superscript also appears somewhere as a left- 
hand index (for instance, in (— l) h a^b\ K c the sum is made over A but 
not over k ) . 

Let a hexahedral cell be given by its eight vertices. Define then edge 
vectors (i/e)„ = o,...,n, node vectors ( M 6)^=0, 1,2, and face vectors ( t /),,=o,....5, 
using the labelling scheme of figure 0a 

» b ■= \ El (4,+,) e /x = 0,1,2 

(12) (-1)' 

1 ' and '/ := —4— ( C 8+20 e +( 9+ 2( l+( -i>'» e) A 

A ((4+20 6 +(5 + 20 e ) L = 0,-..,5, 

with all indices understood cyclic modulo 12 and A denoting the cross 
product in M 3 . 




At every cell face b G {0, 5} and for any given r G M+ the following 
time shifted finite differences of Z in directions M 6 (fj, = 0,1,2) form a 
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vector valued function 
(13) N B Z, t {t) : = 



2(-iy(Z n lt _ T/2 - t Z*\ t ) if /i = [i/2] 
( - 2 »ZP ), t _ T if/i^[i/2] 



( [ x ] denotes the integer part of x G R). The time increments are chosen 
conform with the updating conventions of DSC schemes (as will be seen in a 
moment) and are consistent. In fact, in the first order of the time increment 
r and of the linear cell extension, the vector t V B Z in the centre point 
of face l approximates the scalar products of the node vectors with the 
gradient \7Z . Let, precisely, for a fixed centre point on face i and e G M. + 
the e-scaled cell have edge vectors t e~ : = e L e . Let also ,,V B Z^ denote 
function l|13|) for the e-scaled cell (with node vectors = e^b). Then at 
the fixed point holds 

(14) < „&,grad(Z) > = „b ■ VZ = lim lim Kv B ~Z^, 

as immediately follows from the required C 1 -smoothness of the field Z. 

To recover the gradient VZ from l|13(l in the same order of approxima- 
tion, observe that for every orthonormal basis { v u) v =o,. ..,m-i of M. m or C m , 
and for any basis (^6)^=0,. with coordinate matrix /3£ = < v u , ^b >, 
the scalar products of every vector a with ^6 equal 

(15) < fj), a > = 2^ _ < fj,b, v u > < v u, a > = (3^ a v 





0V = WY 



(at the right-hand side - and henceforth - observe Einstein's convention), 
hence 

(16) a v =Hc£, with (7?) := (G^)*)" 1 • 
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In words: The scalar products of any vector with the basis vectors trans- 
form into the coordinates of that vector with respect to an orthonormal basis 
v u by multiplication with matrix 7 = , where (3^ = < v u, ^b > , 

i.e. /3 is the matrix of the coordinate (column) vectors with respect to 
the given ON-basis v u , and 7 its adjoint inverse. 

This applied to the node vector basis ^b and (|14J) yields the approximate 
gradient of Z at face l 

(17) ,.VZ„ = 7^ t V B Z M . 

The scalar product of the gradient with face vector L f v = < L f, v u > , 
v e {0, 1,2} is thus 

(18) t s = J ■ Nz = j^tf N B z, L = ^ y B z, . 

Continuity of the gradient at cell interfaces yields linear updating equations 
for Z p on the two adjacent faces. In fact, for any two neighbouring cells C, 
X with common face, labelled t in cell £ and n in x > continuity requires 

(19) IS = -IS. 

Substituting l|18|) for JS and *5 and observing the time shifts in 113fl 
provides the updating relations for Z p at the cell interfaces. 
To derive these explicitely, we first introduce the following quantities L z p ' 71 , 
(l = 0,...,5; n = 0,1,2) 



(20) 1*2 (t) 



2(-iy Z n \ t if A* = [i/2] 

(2^+iZ'P - 2 „ZP) lt _ T/2 else 



which in virtue of yields L z p = (p, Z) = ,<r _1 ) = 

= t ;z™ I Z o t <T _1 , where t <r : n 1— > p denotes the nodal shift pertinent to 

face 1 . In particular 

(21) **[!/!!](*) = 2(-l)\^| t , 
which together with (|2U|) for /it 7^ [t/2] is consistent with 

(22) ^(i + r/2) = -l( 2tt+lZ P + 2)iZ P)( t ). 

From (US [H HS1 El) follows that 

»S|t-K = ^(^| t+T/2 - 2(-l)''^/ 2 ] t Z"| t+T ) 

= ^(^| t+ r/2-^ /2I ^l*+-) ■ 

This, with (|19I20|I and the continuity of Z, i.e. = implies 

(24) c , f , _ ^:z»(t + T/2) + ^ ^»(^ + r/2) 
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and for completeness we agree upon setting l zZ(t + T) := L z™(t + r/2), 
for fi 7^ [i/2] . Note that the latter relations contain a slight inconsistency, 
in that continuity might be infringed - which is yet circumvented by taking 
the arithmetic means of the two adjacent values. In fact, our agreement 
doesn't do harm, since any discontinuity disappears with mesh refinement. 

We have thus defined complete recurrence relations for z p ( given z " 
on the former reflection step), which at the same time determine on face l 
the field components and their gradients 

(25) NZ V = tfrf, 

and which essentially constitute the connection step of the algorithm. 
Nodal gradients are similarly, yet more simply, derived using 

v B z;;{t + T -) := { W ZP - 2fl ZV)(t); li = 0,1,2 

in the place of (|13J) and then again l|17JI. With the node and cell-boundary 
values and gradients of T and u , the nodal updating relations for these 
quantities are directly derived from equations ( llOllllI ) in sect. [3]. For equa- 
tions (|10() (without the convective term) this has essentially been carried 
out in |Helj . sect. 5, and the procedure remains straightforward. 

The obtained updating relations are explicit and consistent with near- 
field interaction (only adjacent quantities enter). So, they can optionally 
be transformed into scattering relations for incident and reflected quantities 
© along the guidelines of section [21- with established advantages for the 
stability estimates |He2| . 

5. Pressure 

Pressure is a known subject sui generis in Computational Fluid Dynamics 
|ATP| |MeSt| |CDN| already insofar as pressure fluctuations typically do not 
match the time scales of heat propagation and fluid flow. Pressure fluctua- 
tions are related to acoustic waves, the net effect of which can be important 
and usually has a strong impact on computational stability |LeVeque| . 

Pressure fluctuations play also a key role in the following procedure, 
which is known as divergence cleaning in Magnetohydrodynamics [ibid., 
p. 128] and ensures mass conservation in the present context. 

For Boussinesq-incompressible fluids, conservation of mass simply re- 
quires divergence-free flow, div u = , i.e. in integral form, using Gauss' 
theorem, = div u dV — f 9 * u ■ dF . Since equations (0 00) in sec- 
tion |21 do not a priori ensure this, additional arrangements must be made 
- which is done by means of pressure. 

In a successive overrelaxation (SOR) routine, interposed between the 
connection and reflection steps of the iteration cycle, firstly the (discrete) 
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right-hand side integrals Ig^ = J g( . u ■ dF are computed, and then the 
pressure p which compensates Ig^ so that 



Indeed, taking Z of the preceeding section as the pressure, it follows from 
( 1201 1221 ) that equations l|26[) ( in discrete form with sums over the cell faces, 
of course) yield a unique solution p n = Z n for every cell, given the right- 
hand side integral Iqq . Note that we are actually solving Poisson's equation 
Ap = (qoq/t) divu in integral form. 

With the new face pressure gradient, extracted from ([^OlISlEiD) > the ^ ace 
values of u are updated as u — (t/qoo) gradp. 

After each SOR cycle, continuity of grad p at the cell faces must be 
re-established using (|19|) . i.e. by updating the port values of p according to 
(I23[l . The chain of processes is reiterated until Ig^ < e for a suitable 
bound e (which happens after a few iterations for appropriate choices). 

6. Coaxial line 

To illustrate the approach in a stalwart application, we display the results 
of simulations with coaxial line RL230-100 under high power operating con- 
ditions ( realistically inferred from a ion cyclotron resonance heating ICRH 
experiment in plasma physics). 

The inner and outer conductors of diameters 100 mm and 230 mm are 
made of copper and aluminium, respectively, and the rigid line is filled 
with air at atmospheric pressure. We have simulated the heating process 
from standby to steady state CW operation, at frequency 100 MHz and 160 
kW transmitted power, for horizontal position of the line and with outer 
conductor cooled at 40 degrees Celsius. 

Figure 0] b displays the computed air flow profile (vertical section) in 
steady state, which is attained some minutes after switch-on. Visibly, the 
natural convection pattern is nicely developed. 

The computations have been carried through with a 3D-mesh of 10 layers 
in axial direction, over 200 millimeters of line, the transverse cross section of 
which is displayed in figure 01 a. At the metallic interfaces no-slip boundary 
conditions are implemented and free-slip conditions at all other boundaries. 

Simultaneously, a Maxwell field TLM algorithm was run to provide the 
heat sources. 

7. Conclusions 

A prototypical implementation of the Oberbeck-Boussinesq approxima- 
tion to viscous flow has been presented in this paper, which demonstrates 
the fundamental fitness of the DSC approach for fluid dynamic computa- 
tions. In this respect, at least, (this was recently called into question by a 



(26) 
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TLM expert) DSC schemes significantly transcend the range of application 
of the TLM method, from which they descend. 

A next natural step in the line of this study is, of course, the imple- 
mentation of turbulence models which are compatible with the Boussincsq 
approach, such as the k — e model |ATPj , first of all. - We hope this paper 
stimulates some interest into joint further investigation in that direction. 
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